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In this article, we propose a new method to compute the effective properties of non-linear 
disordered media. We use the fact that the effective constants can be defined through the minimum 
of an energy functional. We express this minimum in terms of a path integral allowing us to use 
many-body techniques. We obtain the perturbation expansion of the effective constants to second 
order in disorder, for any kind of non-linearity. We apply our method to the case of strong non- 
linearities (i.e. D — e{E^)'^^^ E, where e is fluctuating from point to point), and to the case of weak 
non-linearity (i.e. D — eE -\- x{E^)E where e and x fluctuate from point to point). Our results are 
in agreement with previous ones, and could be easily extended to other types of non-linear problems 
in disordered systems. 
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I. INTRODUCTION 

The study of the properties of linear heterogeneous media (such as composites, suspensions) has been the subject 
of an intense activity for already fifty years (see the reviews and j|]). More recently there has been a great interest 
in non-linear media j^]. The non-linearities appear in strong applied fields. In that case, one must expand the local 
generalized susceptibility in powers of the applied field. If the system is invariant under inversion symmetries one 
can write a relation between the displacement field D and the electric field E in the form 

D = eE + xiE^)E (1.1) 

or in the metallic case 

j = aE + x{E^)E (1.2) 
In some cases, the material does not possess a linear regime, and the local equation becomes 

D = x{E'^T''^E (1.3) 

, with K > — 1 (and a similar equation in the metallic case). In the former case, we can speak of 'weak' non-linearity 
^ (WNL) and in the latter one of 'strong' non-linearity (SNL). In the preceding equations, the quantities e and x are 

\ position-dependent and fluctuate randomly from point to point. 

It has been shown Q that the WNL case is connected to conductivity fluctuations, the so-called flicker noise It 
has also been shown Q], that in the limit of weak non-linearity (i.e. xE"^ <C e), the effective non-linear susceptibility 
X* is related to the fourth moment of the local electric field computed in the linear problem (i.e. X = 0). This 
shows in particular that the effective non-linear susceptibility is very sensitive to the micro-structure and that one 
should be careful in building effective medium approximations [|| . This result allowed Stroud and Hui to compute the 
effective susceptibility to first order in impurity concentration. In addition, Zeng et al Q used this relation between 
the fourth moment of the electric field and the effective non-linear susceptibility in order to propose an effective 
medium approximation (EMA). This approximation relies essentially on the 'factorization approximation' (which 
assumes that the local electric field is Gaussian distributed, a questionable hypothesis, especially near the percolation 
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threshold). The critical behavior (near the percolation threshold) was investigated in details by Levy and Bergman 
and numerical simulations were done in In the case where the non-linear term behaves as x(£'^)''/^£', Hui 

and Chung gave a generalization of the effective medium approximation ||T^. 

In the SNL case, the first important result is due to Blumenfeld and Bergman ijllj, [|l2| who derived the second 
order perturbation theory of the effective non-linear susceptibility for three-dimensional composites. This result was 
extended to any dimension d Ponte-Castaiieda proposed a variational approximation in order to establish 
optimal bounds. More recently, EM As (and numerical simulations) were developed in this strongly non- linear case 
II, 0, @. 

For both cases, the EMA reproduces at least qualitatively the numerical results, but there is still a lack of clarity in 
the approximations used. The method we propose here could serve as a starting point for alternative derivations of 
EMA's, and possesses the advantage of not being limited to a special kind of non-linearities. Moreover, we believe 
that this method could be applied to the celebrated random fuse network problem (see for example the review pof). 
The essential point in our method is that one can define the effective constant through the minimum (with some 
constraints) of the stored (or dissipated) energy. We then express this (constrained) minimum in terms of a path 
integral, allowing us to apply usual methods from many-body theory. We derive the second-order perturbation theory 
for any types of non-linearities (leaving more refined approximations for the future [pO[). We illustrate our results on 
both the SNL and WNL cases. 



In the next part, we present the basic equations and derive the second-order perturbation theory. In (IIB) and (11 C), 
we apply our method to the SNL and WNL cases. 

II. CALCULATIONS 

A. Basic equations and second-order perturbation theory 

The effective constants are usually defined through the following relations. In the WNL case, we have 

(D) =e*Eo + x*iE^o)Eo (2.1) 

where Eq is the average field, while in the SNL case 

{D)=x*{Eir^'E, (2.2) 

We first note that there is an alternative definition of the effective constants (see [|l^ and references therein). Let us 
denote by cr^, i = 1, d the generic zero-divergence d-dimensional field. In the dielectric case a is the displacement 
field (Ti = Di whereas in the conducting case it is the current ai — ji ■ The stored energy W (or dissipated energy in 
the conducting case) in the system is given in terms of the energy density w by 

W ^ ^ A'^xw{a{x)) (2.3) 

Hereafter, we take the volume of the sample equal to one. For heterogeneous materials, the energy density depends 
on random quantities fluctuating from point to point. We will specialize our discussion to binary disorder for which 
the parameters can take only two values (our method can easily be generalized to other types of disorder). We thus 
assume that the energy density, at each point, is distributed according to the following distribution probability 

P{w — w{cr{x))) = pS{w — wi{cr{x))) + qS{w — W2{<7{x))) (2-4) 

The alternative of solving Maxwell's equations is to minimize the total energy W subjected to the two constraints 
V • cr = and ai = Ti (the bar denotes the spatial average). As is the case for the free energy in disordered 
systems, one expects this minimum to be self-averaging, allowing us to compute its average over disorder 



(,mm -=T 

Vcr = 



A'^xw{a)) ^W*(T) (2.5) 



where the brackets (•) denotes the average over disorder and W* denotes the stored energy in a homogeneous medium 
characterized by effective constants. The problem thus reduces to the calculation of the average of the minimum of 
a functional of the field subjected to the constraints of zero divergence and fixed mean value. The main point is to 
rewrite this constrained minimum as a path integral 
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min -=T / d'^xw{a{x)) = lim^^oo - / T^'^Si^ ■ a)S{a ~ T)e-''/ (2.6) 

V -=o J P J 

The minimum can thus be interpreted as the ground state energy associated to the partition function Z given by 

(2.7) 



We have to compute the average of the logarithm of (2.7). In order to compute this quantity, we introduce rephcas 



nl and use InZ = Um„^o 



We thus have 



W* hm/5^oolim„^o 



1 (Z'' 



P n 



(2., 



Note that the order of the hmits is important. The rephca method rehes on the fact that one can easily compute 
(Z") for n integer, and then take the limit n 0. The main quantity to study here is thus 



/n 

which for a binary disorder reads 

(Z") = Jf[ Va"SiV ■ a")5(a° - T)exp J d'^a;ln(pe-^S:=i gg-^'SLi -=(-"(-))) 



a=l 



(2.9) 



(2.10) 



Note that in this expression, a cut-off should appear; it wil l only play a role in eq. ( 2.19D , and we will not write it 
in the intermediate expressions. We can now expand ( 2.1C ) to second order in dis order. The fields a"{x) fluctuate 
around their mean value T. To zeroth order, we take a°'{x) = T in equation ( p. 10 ), leading to 



W*iT) ^ (w(T)> 



(2.11) 



In order to set up a perturbation theory to second order, we write a^i^x) = T + e" and expand the exponential in 
(2.10)in powers of e. We obtain (repeated indices are summed) 



/n 

where the tensor Af is given by 



where w[ denotes 4^ at cr = T (and w"- 



*J d(7i da J 



at CT = T). The brackets [A] denote 



[A] 



(A 



— f3nw{T 



^g-/3nto(T)^ 



(2.12) 



(2.13) 



(2.14) 



In the following, we will only need the zeroth order in n of these quantities [A\, which is [A\ = {A) + 0{n). 
At this stage, we have to check that the tensor M is positive definite. Such is the case when the term with the P 
factor is sufficiently small. This term is related to the variance of the disorder, and this requirement coincides with 
the definition of the low contrast limit. We can perform the path integral ( 2.12 ) using the Fourier representation (see 
appendix A for details), and we obtain 



(2.15) 



(where inessential factors were omitted). The first trace is over both space indices and replica indices (a, 7) and 
the second over replica indices only. 

Without loss of generality, we can assume that the energy density is a function of cr^: w — w{(j'^). In this case, the 
tensor M can be written as 
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<5„T,(%A-BT,T,)-/3CT,T, 



(2.16) 



Let us note that, due to spatial isotropy, the only tens ors th at can appear are 5ij and T^Tj. 

One can then compute the determinants in expression (2.15) and one obtains (for n — > and to the first order in the 
variance of disorder) 
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trlnAf ~ (d - l)In(^) + ln(A - BT^ 



/3CT2 
A-BT^ 



(2.17) 



and 



^iAn{kM ifc)~ln(^ A{A - BT^^ ^ [A - BT^){{A - BT^)k'^ - B{k ■ Tf) 
The effective energy is given by the dominant term in /3 



(2.18) 



W*{T) = (w(T)) - 



1 CT2 



2 A - BT2 



1 - 



A 



—y 

>\\d 



{k ■ T)2 



A - ST2 {2AY fc2 + jj^k ■ t) 



(2.19) 



where U = j^^'^yi and we have explicitly written the cut-off A in Fourier space. 

The sum over k can be transformed into an integral (expressed in d-dimensional polar coordinates) 



SdA'' Jo 
Sd-i 

~s7 



dkt 



d-1 



d9,...ded-2 



1 + C/cos^ 



d(?sin 



d-2/ 



COS 



where Sd is the surface of the d-dimensional sphere Sd 

1 CT'^ 
W*{T) = {w{T))--- 



1 + ucos^e 



2 A- BT2 



r(d/2) 
1 



We finally find 

A Sd- 



A-BT^ Sd 



-I{d, U) 



(2.20a) 
(2.20b) 



(2.21) 



where 



I{d, U) 



d6isin''-2g|. 



cos 



1 + Ucos^e 



(2.22) 



These expressions (2.21) and (2.22) are our main results. They represent the expansion of the effective constant to 
second order for any type of non-linearities. 

In the next secti on ( 11 B ) , we apply this result to the case of SNL and we will recover the results of Bergman and Lee 
|0|. In section (II C), we apply our result to the WNL case. 



B. Application to the strongly non-linear case 

In the case of SNL, where D = xix){E^)'^/^ the energy density is given by w{D) — a{x){D^Y^^ where v — ^^ij and 

a — (x) ^ ■ We recall here that in this method we have to express the energy density in terms of the divergence 
free field (i.e. D) and not E. The tensor M here reads 



where 



M;^^ = ,5„^(A% - BT,T,) - pCT^Tj (2.23) 



A = v{a)T'-'^ (2.24a) 
B = v{2 - v){a)T''-'^ (2.24b) 
C ^ v"^ {5a^)T'^''-^ (2.24c) 
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and U = 'jj^ = K. Here and in the following {5XY) = {XY) — {X){Y). We see that M is positive definite, only if 
f3C is small enough, which is the condition for a low disorder expansion. We obtain (using ( ^.21 ) 



W* 



2(^-1) (a) 



2\ r 



V -I bd 



(2.25) 
(2.26) 



Given that a* — (x*) ''+^ , we find the following perturbation expansion for x* for any space dimension (see appendix 
B for details) 

This is the d-dimensional result in the SNL case. It is the same expression as in fl^ ] (this can be seen using the 
change of variables u — cos9, and = 1 — ^ ^, . For d — 2, the result (2.27) reads 



and for c? = 3 



X = X - — ^ — 1 - — ) 
(X) 2k VI + k 



X = ix) ~ 7\ — 1 - ^arcsmJ— -) 

(x) 2k V + 1 



(2.28) 



(2.29) 



C. Application to the weakly non-linear case 

In the WNL case, the energy density is w = e{x){E'^) + x(a;)(£'^)^ which in terms of D gives w{D) = a{x){D'^) + 
h{x){D^Y +0{D^) with a{x) = l/£{x) and h{x) = —x{x)/s'^. We note that although b{x) < 0, there are higher-order 
terms which are positive and guarantee the convergence of the integrals. We note that in this case we can go beyond 
perturbation theory but we will present these results elsewhere |20| . 
In this case, we find 



A =(«)-(&> 
B = 2{b) 



(2.30a) 
(2.30b) 
(2.30c) 



and U — 

In general, the effective energy will not be a polynomial of the form W* — a*T^ + &*T'*. 
effective coefficients, we have to expand the quantity W* to fourth order in T. We obtain 



In order to identify the 



1. 



■(1-3) 



(2.31) 



and 



b*^{b)-2{l--) 



(3- 



d+2 



{b){6a^) JSab). 



(2.32) 



We see in these expressions that for d = 1, we obtain a* = (a) and 6* = (6). This is the exact result since for d = 1, 
the constraints V • ct = and a = T imply a = const — T. 

Knowing that a = 1/e and = x/^"'j we obtain after straightforward but tedious calculations 



£ - 



1 



(2.33) 
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and 



X* - (x> 



4(2-i)iM 



2( 



10rf2 + l5rf- 16 (fe 



d{d + 2) 



(2.34) 



These are the second-order perturbation results for the effective coefficients in the WNL case. We note the surprising 
resuh that, at this order, e* does not depend on moments of x ^nd that x* does not depend on {Sy^). This suggests 
that the effective medium result for e* should be independent of the non-linearity (in this weak non-linearity limit) 
and that x* is not determined self-consistently but is a function of e* (this seems to be the case in a self-consistent 
resummation of the perturbation theory pO|). 



III. CONCLUSION 



In this article, we proposed a new method to compute effective properties of disordered non-linear media for any 
type of non-linearities. In the case of strong non-linearity, we recover the known perturbation results and we present 
the second order perturbation result for the weakly non-linear case. We would like to emphasize that our method is 
general and could be used in many other situations (such as plasticity of porous media or the random fuse network) . 
Moreover, this framework could be used as a starting point for more refined approximations. 
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APPENDIX A: 



In this appendix, we compute the path integral in expression (2.12). For this purpose, we will use the Fourier 
representation of the field ta{x) (we keep the same notation for the field and its Fourier transform) 



(Al) 



The integral ( 2.12 ) can thus be rewritten as 
d'^?„(0)e-^^^(")-*^-^^(") 



By exponentiating the (5-functions, one finds after integration the result (Eq. (2. IE)) of the main text 



7n\ ^, gln(e 



-^„™(T)< -iV trlnM -iV trln(fcM"^fe) 



(A2) 



(A3) 



APPENDIX B: 



We explicit here the relation between the perturbation expansion of a and that of x- We look for an expansion of 
X* of the form 

X* = (x)-^.9x (Bl) 
We expand this expression in power of the contrast A = X2 ^ Xij and find 

X* XI + gA - P^g^ (B2) 
Xi 
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Now, assuming that the perturbation expansion for a is of the form 



a = (a) 



-9a 



(B3) 



this leads to (to second order in A) 



1/(k+1) 



Xi(« + 1) + X?(« + 1)^^ 2 ^^"^ 



(B4) 



On the other hand, the relation between the functions ga is g^^ is found using a* = (%*) Expanding this 

relation to second order in A, we find 



a ~xi 



qA 



gA2 ^-?«: + 2, 



Comparing expressions (B4) and (B5), we obtain 



9x^o 



1 K + 2 ga 



2 K+1 K+1 



(B5) 



(B6) 



This relation allows to relate g^ and ga- 
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